from astropy.io import fits
import numpy as np

def gaussian(r,sigma):
    return 10000*np.exp(-r**2/2/sigma**2)
im = np.zeros((4000,4000))
xx,yy = np.meshgrid(np.arange(4000),np.arange(4000))
xx -= 2000
yy -= 2000
rr = np.sqrt(0.75*xx**2+0.75*yy**2-0.5*xx*yy)
im = gaussian(rr,200)
print(im)
im1 = np.zeros((40,40))
for i in range(40):
    for j in range(40):
        im1[i,j] = np.mean(im[100*i:100*(i+1),100*j:100*(j+1)])+2
fits.writeto("testim.fits",data=im1,overwrite=True)
